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Abstract: Determination of pressure distribution on airfoil 
profiles enables to understand the dynamic flow characteristics 
around the surface of wings or blades. Computational methods 
are widely applied in the field of aerodynamics to predict the 
aerodynamic characteristics such as the pressure , lift and drag 
coefficients on an airfoil. In the present paper, the pressure 
distribution of NACA airfoil profiles is calculated using the 
numerical panel method for 2D lifting air flow conditions. The 
analysis of airfoil geometry subjected to various AO A (angle of 
attack ) including the stalling angles is considered to observe 
the variation in pressure coefficients along the chord. The zero 
lift AO A for the profiles are also evaluated in order to assess 
influence of thickness to chord ratios on airfoil characteristics. 


line or camber of the profiles. The normalized coordinates are 
obtained by the dividing the x and y values with the chord of 
airfoil which can be then used for analysis. The profiles of 
NACA 4412, NACA 16-006, NACA 0024 and NACA 66-018 
are shown below. The individual profiles can be differentiated 
using the distinct parameters as leading edge radius, trailing edge 
angle. The airfoils have finite thickness at trailing edges which 
affect the aerodynamic properties in terms of noise and stalling 
behavior. The leading edge radius and trailing edge angle for 
NACA MPXX 4 digit airfoils are [1, 3] 

t = L 1019 ' [;f w 

5 te = 2 tan -1 [1.16925 -] ...(B) 
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1. Introduction 

NACA airfoil profiles are widely used in the aircraft industry for 
producing lift on the wing span. The geometry of airfoil dictates 
the airfoil performance and determines its relevance for a 
specific application. The present analysis deals with the analysis 
of pressure distribution over the airfoil surfaces along with the 
lift and drag characteristics. They are commonly expressed in 
the form of drag polar for the NACA series and represent the 
aerodynamic performance at various flow field conditions. The 
Reynolds number is often used to characterize the flow 
conditions and predict the behavior of airfoils at various 
operating flow conditions both for the viscous and in viscid 
fluids. Typically the results from the experiments serve as the 
emulation tool in the aircraft industry for comparison with actual 
performance of the wing. Panel methods are modern numerical 
techniques which are quick to execute and predict accurate 
results compared to the traditional experimental methods that are 
cumbersome procedures and time consuming. 


The thickness and mean line distributions for the 4 digit airfoils 
are expressed in terms of the t/c ratio, and slope of camber line is 
obtained by differentiating the camber line function, y c [6, 1, 7] 
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Where the a 0 , ai a 2 , a 3 & a 4 are the coefficients describe the shape 
of airfoil and its distribution of ordinates along the chord line. 


2. Airfoil geometry 

NACA airfoils are designed during the period from 1929 and 
1947 by Eastman Jacobs at NASA Langley field laboratory [2, 
8]. The airfoil geometry for most of the NACA profiles can be 
divided into x-coordinates known along th Qchord line and y - 
coordinates known as the ordinates, which represent the mean 
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Figure 1 : Geometry of NACA 4 digit, modified & 6 series profiles 


3. Computational panel method 

Traditional methods for modeling flow around slender bodies of 
any shape include potential flow which utilizes the superposition 
of source and sink on x axis and in uniform distributed flow. 
However, the theory does not predict accurate values for flow 
whose leading edge has rounded shapes. Basic panel methods 
were developed by Hess and Smith at Douglas aircraft in late 
1950s [2] for aircraft industry. Panel methods model the 
potential flow by distributing sources over the body surface. A 
source is point at which the fluid appears in the field at uniform 
rate while a sink is point which disappears at uniform rate, m 3 /s. 
Each source or sink has specific strength denoted by circulation, 
r. Simple 2D uniform lifting flows [1] can be described using 
the following equations 


cp = by — ax ... stream line function 
0 = bx + ay . . . Velocity potential function 


The resultant velocity V, at any point along the flow direction 
can be written as 


Via 2 + b 2 1 


• Panels are defined with unit normal and tangential 
vectors, n , f . 

• Velocity vector, denoted by v^are estimated by 
considering the two panels, i & j the source on the panel j which 
induce a velocity on panel i. The perpendicular and tangential 
velocity components to the surface at the point I, are given by 
scalar products of . n and . t 

• The above quantities represent the source strength on 
panel j and expressed mathematically as 

Vij.n = tfjNjj 
Vjj-t = GjTjj 

Where INE and T^are the perpendicular and tangential velocities 
induced at the collocation panel i and known as normal and 
tangential influence coefficients. The surfaces represented by the 
panels are solid and the following conditions are applied for the 
normal and tangential velocities at each of collocation points 
consisting of sources strengths, vortices, and oncoming velocity, 
U. 
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Sources Vortices Oncoming flow 


The above system of linear algebraic equations are solved for the 
N unknown source strengths, a i9 using matrix system and 
expressed as 


The point source or sink is distributed uniformly in all 

directions of the flow field and obey the continuity equation and 
irrotational motion everywhere except at the point itself. The 
total velocity potential for numerical panel method can be 
rewritten as follows [2] 

0 = Ux + g 

Where, p is the perturbation potential away from free stream 
conditions. The airfoil geometry is discretized into finite number 
of panels over the surface. The panels are represented by the 2D 
shape of the surface by series of straight line segment [2, 4]. 
The following procedure describes the calculation for 2D lifting 
flows. 

• Numbering of end points or nodes of the panels from 
1...N 

• The center points of each panel are chosen as 
collocation points. The boundary condition of zero flow 
orthogonal to surface is applied to the points. 


M.a = b ... (4) 


Where N is an N+l x N+l matrix containing the N^and is 
column matrix of N elements and A is the column matrix of N 
elements of unit normal velocity vectors. Matrix inversion 
procedures available in MATLAB are applied to solve for the 
source strengths using the above system of equations. The 
pressure acting at collocation point i is given by the Bernoulli 
equation as [2, 4, 5] 


r = 1 — 

'-'pi ± 


v Ti 

U 


...(5) 


Where v Ti the tangential velocity vector is determined using the 
influence coefficients. The influence coefficients are important 
for panel method in order to determine the pressure distribution 
over the surface of the any given airfoil coordinates. The airfoil 
trailing edge presents a unique condition for the flow field 
parameters. Using panel method, the following criteria is applied 
for the stream lines around the airfoil 
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• The streamlines leave the trailing edge with a direction 
along the bisector of the trailing edge angle. 

• The velocity magnitudes on the upper and lower 
surfaces near the trailing edge of airfoil approach the same 
limiting values. 

• The trailing edge angle is modeled as the stagnation 
point for finite value of trailing edge angle hence the source 
strength must be zero at the trailing edge. 


The above assumptions are known as the Kutta condition which 
is essential for the successful evaluation of velocity vectors and 
pressure for 2D uniform flows. It can be written in the algebraic 
form of equation as 

E°j T ,,j+Y T ,,N+i + U-t i = 


j=l 


f n _ ^ 'N 

Vj =1 J 



It must be noted that the velocity components induced at any 
point P due to sources on panel centered at point Q can be 
expressed mathematically as scalar product of vector form 


v P Q = v xQ .tj + VyQ . fij ...(7) 
Njj = VpQ.fi; ...( 8 ) 

Tjj = VpQ.tj ...(9) 


Therefore, the resulting velocity components along with known 
source strengths and influence coefficients are added for each 
panel in order to obtain pressure distribution over the airfoil 
surface. The number of panels used (order) in the simulation 
ranged from 50 - 70 for airfoils in the MATLAB routine foil.m. 
Depicted below is the flow chart for the computational panel 
method described in previous section. 

Flowchart: Computer routine for panel method 


NACA 0015 Airfoil pressure distribution 



Figure 3: Pressure distribution of NACA 0015 airfoil 

The following are the list of NACA airfoils for which 
the panel method is used to calculate the pressure distribution. 
The ideal or design lift coefficient is calculated for the zero angle 
of attack using the formula [1,3] 


Cnh 


OCi = 


...(9) 


2tu(1 + a) 

Where Cli is design lift coefficient, h is the function that 
describes the uniform load along the mean line and the type of 
mean line used in the airfoil. The thickness and camber is written 
mathematically as 
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Figure 4: Angle of attack at zero lift along the wing span 



4. Results and Discussion 


-.100 = % cambe; 100 = % thickness ... (10) 
c c 

Fig. 2 shows the computed values for pressure distribution of 
NACA 0015 airfoil. It has 15 % thickness and 0 % camber and 
identical to the symmetrical airfoil. The 4 digit and 4 digit 
modified airfoils have same thickness distribution below and 
above the mean line; however they are different in the mean line 
or camber defined by the slope at the leading edge. It can be 
noted that the symmetric airfoils have good stall behavior and 
tend to exhibit high pitching moment coefficients. They also 
present a high drag and insensitive to roughness. Some of the 
common applications where they are used include helicopter 
blades, missiles and shrouded turbines in general aviation 
industry. 
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Table 1: NACA Airfoil series nomenclature 


SI. 

Mn, 

NACA 

\ iffnil 

Thickness 

cl 

Maxim 


1 

0015 

15 



2 

0024 

24 

- 

- 

3 

4412 

12 

0.04, 0.4 

- 

4 

66018 

18 

0.06, 0.6 

- 

5 

16006 

6 

0.06, 0.6 

- 

6 

64206 

6 

0.04, 0.4 

0.2 


Fig 2 & Fig 3 shows the NACA OOXX series profiles which 
produce different characteristics at range of angle of attack i.e. -1 
deg to 3 deg respectively. The upper (suction) and lower 
(pressure) surfaces are both closer in case of NACA 0015 while 
NACA lower surface pressure coefficients are lower compared 
to upper surface in NACA 0024 as the angle of attack is 
changed. Further, with the change in onset flow velocity, i.e with 
increasing Reynolds number, the relative values for pressures on 
the surfaces show significant variances. Towards the trailing 
edge, the pressure gradient suddenly drops where the velocities 
at the upper & lower surfaces reach nearly the same limiting 
values per the kutta condition. 


NACA 0024 airfoil ■ pressure distribution 
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Figure 5: Pressure distribution of NACA 0024 airfoil 

From fig 4 the NACA 4412 profile shows the significant high 
drag and low maximum lift coefficient and as the incidence 
angle is changed the upper and lower surface pressure values gap 
narrow. Further, at higher Reynolds number, or negative AoA 
[9] it shows the pressure reversal of upper (suction) and lower 
(pressure) surfaces can be observed. Fig 9 shows the comparison 
of pressure distribution indicate that the NACA 6 series exhibit 
high maximum lift coefficients compared with the 4 digit and 
modified 4 digit series. 


NACA 4412 Airfoil pressure distribution 



Figure 6: Pressure distribution of NACA 4412 airfoil 


NACA 2412 Airfoil 



Figure7: Pressure distribution of NACA 2412 airfoil 

NACA 16-006 airfoil - pressure distribution 
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Figure 8: Pressure distribution of NACA 16-006 airfoil 
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NACA66- 018 airfoil pressure distribution 



Figure 9: Pressure distribution of NACA 66-018 airfoil 


NACA 64206 airfoil pressure distribution 



Figure 10: Pressure distribution of NACA 64-206 airfoil 
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Figure 11 Comparison of pressure distribution for various NACA profiles 

using computational method. 


The lift and drag characteristics were also computed using panel 
method by resolving the resulting force which is inclined at 
small angle with the AoA on the suction side of airfoil. The lift 
and drag can be expressed as 


C L = C z .cosa— C x . sin a ... (10) 
C D = C z .cosa+ C x . sin a ... (11) 


Where, C z and C x are the coefficients along the Y and X 
directions of the airfoil coordinate system. The values are 
obtained by pressure gradient coefficients of the upper and lower 
surfaces of airfoil and multiplying with the length of panel along 
the Y and X directions for each panel. The drag polar are plotted 
in fig 12 and 13 for the NACA 00XX airfoils. A drag polar is 
any mathematical function which relates the drag coefficient 
with some function of lift coefficient. Most often the drag 
coefficient is dependent upon the Reynolds number and Mach 
number and also expressed in terms of aspect ratio of wing. 
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Figure 22: Lift and Drag Characteristics of NACA 0024 airfoil 



Figure 33: Comparison of CL/CD ratios of NACA OOXX series 
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5. Conclusions 

NACA airfoil study on pressure distribution, lift and drag 
attributes was performed using the computational panel method 
for 2D lifting flows. From the results this method is more useful 
when the computational requirements are low and with reduced 
time to validation of the results. It also predicts the aerodynamic 
characteristics with reasonable accuracy. The number of panels 
defined for the airfoil will affect the accuracy of the predicted 
values and its flow characteristics parameters. It overcomes the 
limitations present in the traditional methods developed for the 
flows over the non-lifting bodies. 
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